#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <math.h>
#include <gsl/gsl_integration.h> 
#include "astro/utilities.h"
#include "astro/constants.h"
#include "astro/cdmModel.h"
#include "astro/cosmology.h"
#include "astro/cdmLens.h"

const int n=128,ni=512;
const double tiny=1e-6;
gsl_integration_workspace * work= gsl_integration_workspace_alloc (n);
const size_t limit = n;
std:: vector<double> dci(n),lzi(n),zi(n);
std:: vector<double> vii(ni),lmii(ni);
const double mmin=1e3,mmax=1e20;

/*
  Initialise vectors for MAH
*/
void iniTables(astro:: cdmModel *cdm,double z0){
  double lz0 = log10(1+z0);
  astro:: fill_linear(lzi,n,0.+lz0,1.+lz0);
  for(int i=0;i<n;i++){
    zi[i]=-1+pow(10,lzi[i]);
    dci[i] = cdm->deltaC(zi[i])/
    	(cdm->growthFactor(1./(1+zi[i]))/cdm->growthFactor(1.)/(1.+zi[i]));
  }
  astro:: fill_linear(lmii,ni,log10(mmin),log10(mmax));
  for(int i=0;i<ni;i++){
    vii[i] = cdm->variance(pow(10.,lmii[i]));
  }
}


/*
  estimate the mass from the variance
*/
double massFromVariance( double v )
{
  if( v>vii[0] )
    return mmin;
  if( v<vii[ni-1] )
    return mmax;
  int i = astro:: locate (vii, v);
  i = std::min (std::max (i,0), int (ni)-2);
  double f=(v-vii[i])/(vii[i+1]-vii[i]);
  return pow(10.,f*lmii[i+1]+(1-f)*lmii[i]);
}


/*
  estimate redshift from \delta_c
*/
double getZfromdc(double dc){
  if(dc>dci[n-1])
    return zi[n-1];
  if(dc<dci[0])
    return zi[0];
  int i = astro:: locate (dci,dc);
  i = std::min (std::max (i,0), int (n)-2);
  double f=(dc-dci[i])/(dci[i+1]-dci[i]);
  return -1 + pow(10,(f*lzi[i+1]+(1-f)*lzi[i]));
}
/*
  mean formation redshift definition by half mass of the main halo progenitor
  extended to the ellipsoidal collapse
*/
double getZf(astro:: cdmModel *cdm, double m0,double z0){
  // from the ellipsoidal collapse model Giocoli et al. 2007
  double q = 0.707; 
  double s0 = cdm->variance(m0);
  double s0_2 = cdm->variance(m0/2.);
  double dc0 = cdm->deltaC(z0)/
    	(cdm->growthFactor(1./(1+z0))/cdm->growthFactor(1.)/(1.+z0));
  double df = dc0 + 0.974/sqrt(q)*sqrt(s0_2-s0);
  return getZfromdc(df);
}
/*
  begin: Routines for the mass accretion history by Zhao et al. 2009
*/
double DlnSigmaDlndeltaC(astro:: cdmModel *cdm,double m0,double z0,double mp,double zp){
  double dlsdlm0 = cdm->DlnSigmaDlnM(m0);
  double S0 = sqrt(cdm->variance(m0))*pow(10.,dlsdlm0);
  double Dp0 = cdm->growthFactor(1./(1+z0))/cdm->growthFactor(1.)/(1.+z0);
  double w0 = (cdm->deltaC(z0)/Dp0)/S0;
  double p0 = 0.5*w0/(1+pow(w0/4.,6.));
  double Dpp = cdm->growthFactor(1./(1+zp))/cdm->growthFactor(1.)/(1.+zp);
  double dcp = cdm->deltaC(zp)/Dpp;
  double incz = (log(dcp) - log(cdm->deltaC(z0)/Dp0))/0.272/w0;
  double pz;
  if(incz<1) pz = p0*(1-incz);  
  else pz = 0;
  double dlsdlmp = cdm->DlnSigmaDlnM(mp);
  double Sp = sqrt(cdm->variance(mp))*pow(10.,dlsdlmp);
  double wp = dcp/Sp; 
  return (wp - pz)/5.85;
}
/*
  mass accretion history by Zhao et al. 2009
*/
double mainProg(astro:: cdmModel *cdm,double m0,double z0,double mp,double zp,double zi){
  double dlnsdlnM = DlnSigmaDlndeltaC(cdm,m0,z0,mp,zp);
  double Dpp = cdm->growthFactor(1./(1+zp))/cdm->growthFactor(1.)/(1.+zp);
  double dcp = cdm->deltaC(zp)/Dpp;
  double Dpi = cdm->growthFactor(1./(1+zi))/cdm->growthFactor(1.)/(1.+zi);
  double dci = cdm->deltaC(zi)/Dpi;
  double DlnS = dlnsdlnM*(log10(dci) - log10(dcp));
  double sigp = sqrt(cdm->variance(mp));
  double sigi = pow(10.,log10(sigp) + DlnS);    
  return massFromVariance(sigi*sigi);
}
/*
  return redshift when a given fraction, mf,  of the main halo progetor
  has been assembled
*/
double getZfram0(astro:: cdmModel *cdm,double m0,double z0,double mf){
  std:: vector<double> mi(n),lmi(n);
  double mp=m0;
  double zp=z0;
  for(int i=0;i<n;i++){
    if(i==0) mi[i] = m0;
    else mi[i] = mainProg(cdm,m0,z0,mp,zp,zi[i]);            
    if(mp<mi[i]){
      std:: cout << " WARNING: mass increases with z " << std:: endl;
      // std:: cout << " I will STOP here!!! " << std:: endl;
      std:: cout << i << "  " << mp << "   " << mi[i] << "  " <<z0 << "  " << zi[i] << std:: endl;
      // exit(1);
    }
    mp=mi[i];
    zp=zi[i];
    lmi[i] = log10(mi[i]/m0);
  } 
  // interpolate to give back the redshift
  double lmf = log10(mf);
  if( lmf<lmi[n-1] ) return zi[n-1];
  int i = astro:: locate (lmi, lmf);
  i = std::min (std::max (i,0), int (n)-2);
  double f=(lmf-lmi[i])/(lmi[i+1]-lmi[i]);
  return -1+pow(10,(f*lzi[i+1]+(1-f)*lzi[i]));
}

/*
  estimate the concentration as function of the time
  when the host halo assmble 4% of the final mass
 */
double getC(astro:: cdmModel *cdm,double m0,double z0){
  double t0 = cdm->time(z0);
  double z004 = getZfram0(cdm,m0,z0,0.04);
  double t004 = cdm->time(z004);  
  return 4*pow(1.+pow(t0/3.75/t004,8.4),1./8.);
} 

int main(){ 
  std:: string filein,str,fileout;
  std:: ifstream fin;
  std:: ofstream fout;
  filein = "INPUT";
  fin.open(filein.c_str());
  if(!fin){
    std:: cout << "INPUT file does not exist" << std:: endl;
    std:: cout << " I will STOP here!!!" << std:: endl;
    exit(1);
  }
  double om0,oml,h0,sig8,z0,w,lminmass,lmaxmass;
  fin >> str;
  fin >> om0;
  fin >> str;
  fin >> oml;
  fin >> str;
  fin >> h0;
  fin >> str;
  fin >> sig8;
  fin >> str;
  fin >> w;
  fin >> str;
  fin >> z0;
  fin >> str;
  fin >> lminmass;
  fin >> str;
  fin >> lmaxmass;
  fin >> str;
  fin >> fileout;
  fin.close();  
  // read in put file end initialise the cosmological model and vectors
  astro:: cosmology co(om0,oml,h0,w);
  astro:: cdmModel cdm(co,8.,sig8);  
  std:: cout << " cosmology: " << std:: endl;
  std:: cout << "     om0="<<om0<<"   oml="<<oml<< std::endl;
  std:: cout << "     h0="<<h0<<"   sig8="<<sig8<< "   w="<<w<< std::endl;
  std:: cout << "    mstar("<<z0<<") = " << cdm.nonlinMass(z0) << std:: endl;
  std:: cout << "   " << std:: endl;
  std:: cout << "   mc-relation will be computed for haloes at z="<<z0 << std:: endl;
  std:: cout << "   " << std:: endl;
  int nh=128;
  std:: vector<double> lm0;
  // double lminmass = log10(20.*3.929*pow(10.,5)*h0);
  astro:: fill_linear(lm0,nh,lminmass,lmaxmass);
  std:: vector<double> m0(nh);
  iniTables(&cdm,z0);  
  fout.open(fileout.c_str());
  for(int i=0;i<nh;i++){
    m0[i] = pow(10.,lm0[i]);
    fout << m0[i] << "   " << getC(&cdm,m0[i],z0) << std:: endl;
  }
  fout.close();
  std:: cout << " output file of the mass concentration relation is in the file:" << std:: endl;
  std:: cout << fileout << std:: endl;
  std:: cout << "   " << std:: endl;
}
